Multilayer graphenes with mixed stacking structure 
— interplay of Bernal and rhombohedral stacking 
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We study the electronic structure of multilayer graphenes with a mixture of Bernal and rhom- 
bohedral stacking and propose a general scheme to understand the electronic band structure of an 
arbitrary configuration. The system can be viewed as a series of finite Bernal graphite sections 
connected by stacking faults. We find that the low-energy eigenstates are mostly localized in each 
Bernal section, and, thus, the whole spectrum is well approximated by a collection of the spectra 
of independent sections. The energy spectrum is categorized into linear, quadratic and cubic bands 
corresponding to specific eigenstates of Bernal sections. The ensemble-averaged spectrum exhibits 
a number of characteristic discrete structures originating from finite Bernal sections or their com- 
binations likely to appear in a random configuration. In the low-energy region, in particular, the 
spectrum is dominated by frequently-appearing linear bands and quadratic bands with special band 
velocities or curvatures. In the higher energy region, band edges frequently appear at some particular 
energies, giving optical absorption edges at corresponding characteristic photon frequencies. 

PACS numbers: 73.22.Pr 81.05.ue,73.43.Cd. 



I. INTRODUCTION 

Graphene multilayers exhibit a wide variety of stack- 
ing arrangements allowed by weak van der Waals inter- 
layer coupling, offering various types of quasiparticles in 
the low-energy electronic spectrum. In three-dimensional 
bulk graphite, there are two distinct crystal configura- 
tions called Bernal (ABAB • • • ) [H-0] , and rhombohedral 
(ABC ABC ■■■) stacking @|| as illustrated in Fig. [Q A 
sequence such as ABC ■ ■ ■ represents the lattice point on 
every layer along a perpendicular axis, where A and B 
are inequivalent sublattices of hexagonal lattice, and C 
is the center of the hexagon. Recently, several exper- 
imental techniques, such as optical absorption [10l - [l3l |. 
Raman spectroscopy [l4], [H| and transmission electron 
microscopy [l6j], have been applied to identify the num- 
ber of layers and the stacking order of graphene multi- 
layers. Few-layer graphene samples exfoliated from bulk 
graphite usually exhibit Bernal structure which is sup- 
posed to be the most stable, but often also display rhom- 
bohedral structure in part [15, 

The electronic band structure of graphene multilayer 
depends sensitively on its stacking structure |15l - l23| . 
In Bernal stacked multilayer, the spectrum consists of 
quadratic bands analogous to bilayer graphene and a sin- 
gle linear band like monolayer [24l - |29| . In contrast, a 
rhombohedral-stacked multilayer has a totally different 
spectrum with a pair of fiat low-energy bands which dis- 
perse as p N with momentum p and the number of layers 

In general, graphitic structures are expected to take a 
Bernal-rhombohedral mixed form as illustrated in Fig. [5] 
The energy spectrum of mixed multilayer was studied for 
some specific few-layer cases [IE IIS III], but general rules 
predicting the electronic properties of arbitrary struc- 
ture are not well known. A single rhombohedral stacking 
fault appearing in Bernal graphite was studied theoret- 



ically [3J], and it supports cubic bands associated with 
the localized states bound to the stacking fault. It was 
also shown [32j that the low-energy spectrum of Bernal- 
rhombohedral mixed multilayer consists of energy bands 
which disperse as p J , and the sum of J coincides with the 
number of layers. This predicts the number of the bands 
belonging to each J, but to obtain the quantitative dis- 
persions and wavefunctions one would need to actually 
calculate eigen energies for every single configuration. 

In this paper, we study the electronic structures of 
general Bernal-rhombohedral mixed graphene multilay- 
ers. We begin by proposing a general scheme in which to 
understand the band property of any given configuration 
without resorting to diagonalizing the full Hamiltonian. 
We view the system as a series of finite-layered Bernal 
sections connected by the rhombohedral-type stacking 
fault, as depicted in Fig. Efa), and treat the coupling 
between neighboring sections as a perturbation. We find 
that the eigenstates near the Dirac point are mostly lo- 
calized in each Bernal section, and the states are approx- 
imated well by those of incomplete Bernal graphites as 
illustrated in Fig. [3Jb) and (c). The energy spectrum 
is then categorized into linear, quadratic and cubic (or 
higher order) bands within the basis of incomplete Bernal 
graphite sections. We also specify several limited situa- 
tions where the states of neighboring Bernal sections are 
strongly hybridized. 

To model realistic experimental systems with many 
layers and, perhaps, a number of stacking faults, we 
develop a statistical approach. To study the electronic 
properties averaged over different stacking configura- 
tions, we analyze the statistics of the velocity of the lin- 
ear bands and the effective mass of the quadratic bands 
which dominate the low-energy spectrum. We find that 
there are some particular frequently-appearing values in 
the velocity /mass distribution, corresponding to finite- 
layered Bernal sections and their particular combina- 
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FIG. 1: Lattice structures of (a) Bernal graphite and (b) 
rhombohedral graphite. In each panel, the right figure is a 
top-view, the middle is a schematic diagram of the lattice 
structure. 



tions. We also compute the averaged optical absorp- 
tion spectrum and find that the absorption edges emerge 
at particular frequencies, corresponding to frequently- 
appearing band structures. 

The paper is organized as follows. We formulate the 
Hamiltonian of mixed multilayer graphene in Sec. II. We 
describe the eigenstates and the energy spectrum of iso- 
lated incomplete Bernal graphite section, and study the 
inter-section mixing effect in Sec. III. We present the 
ensemble-averaged distribution of low-energy band ve- 
locity and effective mass, and the optical absorption in 
Sec. IV. 



II. FORMULATION AND EFFECTIVE-MASS 
HAMILTONIAN 

We consider an V to t-layer graphene system which has 
a unit cell containing Aj and Bj sublattices on the j-th 
layer. Coupling between the j-th and j + 1-th layers is de- 
scribed as either AB or BA stacking, where AB stacking 
is defined as the arrangement in which sites Aj and -Bj+i 
are connected by the vertical interlayer coupling, whereas 
BA stacking as the one in which Bj and Aj + \ are con- 
nected, as illustrated Fig. [5] The entire system is speci- 
fied by a set of indices aj = ± for j = 1, 2, ■ • • , iVtot - 1 de- 
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FIG. 2: Example of Bernal-rhombohedral mixed graphene 
multilayers expressed by (5,7,4,1,3). Intralayer coupling by 
parameter 70 is shown as a vertical solid line, interlayer cou- 
pling by 71 is shown as a horizontal solid line. 



scribing the interlayer connection between j-th and j + 1- 
th layers, where + and — represent AB and BA stack- 
ing, respectively. Bernal-stacked graphene is then ex- 
pressed as an alternating sequence like (+,—,+,—,•••), 
while rhombohcdral-stacked graphene is (+, +,+,•••) or 

(---■••)• 

In a general configuration, an alternative sequence 
(••• ,+,—,+,—,•••) is regarded as a section of contin- 
uous Bernal structure, and a position at which the same 
sign (• • • ,+,+,•••) or (• • • ,—,—,■••) occurs consecu- 
tively is regarded as a rhombohedral-type stacking fault 
separating different Bernal section. The sequence {oi} is 
alternatively expressed as a set of integers 



(N U N 2 ,N 3 ,--- ,N M ), 



(1) 



where 2Vj is the length of the i-th Bernal section, and the 
stacking fault exists between iVj and iVj+i. For exam- 
ple, the sequence (+, — , +, — , — , +, — , — , +) is written as 
(4,3,2). M is the number of separated sections in the 
whole system. The total number of layers in the system 
is given by 



N tnt = 1 



M 



N, 



(2) 



A pure Bernal-stacked multilayer graphene is represented 
by a single number (JV to t — l)j and a pure rhombohedral 
multilayer is by (1, 1, 1, • • • ). 

To describe the electronic properties, we use an 

effective-mass model 0, H, [3oT - l3a | with the Slonczewski- 
Weiss-McClure parameterization of graphite [6]. As the 
simplest approximation, we include parameter 70 describ- 
ing the nearest neighbor coupling within each layer, and 
71 for the coupling of the interlayer vertical bonds. The 
band parameters were experimentally estimated in the 
bulk ABA graphite, for example [6] as 70 = 3.16 eV and 
71 = 0.39 eV. The low energy spectrum is given by states 
in the vicinity of the K% point at the corner of the Bril- 
louin zone, where £ = ±1 is the valley index. If \Aj) and 
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\Bj) are Bloch functions at the K% point, corresponding 
to the A and B sublattices of layer j, respectively, then, 
in the basis of \A 2 ), \B 2 ), ■ • • , the Hamiltonian 

in the vicinity of the valley is 



(a) 



/H Vx 
Vl H V 2 



H 



\ 



Vi Ho V 3 
Ho 



V 4 



(3) 



with 




= -), 



(4) 



(5) 



Here, the in-plane momentum operator is tt = £p x + ip y , 
and p = (p x ,Py) = —ihV. The diagonal blocks, Eq. @, 
describe nearest-neighbor intralayer hopping, and V, 
Eq. ([5]), describes nearest- neighbor layer hopping, v 
is the band velocity of monolayer graphene given by 
v = \/3a7o/2?i, where a 0.246 nm is the lattice con- 
stant of honeycomb lattice. We neglect other hopping 
parameters in the following arguments for simplicity. We 
expect that the neglected parameters introduce relatively 
small corrections such as trigonal warping and electron- 
hole asymmetry, as in multilayer graphenes with pur e 
Bernal H, [H, [23, [3!| and rhombohedral stacking [9l.l33j| . 



III. ELECTRONIC STRUCTURE AND BAND 
CLASSIFICATION 

Unlike a pure Bernal multilayer, the Hamiltonian of 
a general graphene stack cannot be simply decomposed 
into smaller subsystems by a unitary transformation. At 
small momentum, however, we can show that every eigen- 
state is almost well localized in a single Bernal-graphite 
section between rhombohedral stacking faults, so that the 
system can be treated approximately as a set of indepen- 
dent Bernal sections. 

In order to develop an analytical description, we be- 
gin by considering an imaginary system in which the in- 
tralayer hopping is switched off only on the stacking-fault 
layers, as illustrated in Fig. G^a). The entire system is 
then broken into a set of incomplete Bernal graphite sec- 
tions with two inequivalent sublattice sites per layer but 
with one sublattice missing at the stacking fault layers, as 
illustrated in Fig. [3]Jb) and (c) for a middle section and 
an end section, respectively. In the following, we will 
show that at small momentum p -C Ji/v, the spectrum 
is approximately that of a collection of the energy bands 
of isolated incomplete Bernal graphites, except for some 
specific occasions where the coupling between different 
sections is significant. 




FIG. 3: (a) Decomposition of Bernal-rhombohedral mixed 
multilayers into incomplete Bernal sections. Dashed lines are 
connections between neighboring sections, which are to be 
switched off in the decomposition, (b) Decomposed incom- 
plete Bernal section in the middle and (c) at one end of the 
whole system. 



A. Spectrum of incomplete graphite 

We first consider the electronic structure of an iso- 
lated incomplete Bernal graphite of the middle section 
type, as in Fig- EJ^b)- For convenience, we divide all the 
atomic sites into two groups, 'chained' sites and 'free' 
sites, where a chained site refers to a site connected to 
neighboring layers by 71 , and a free site to one which has 
only intralayer connection 70. This incomplete Bernal 
middle section has N + 1 layers that may be numbered 
from j = to N. In this respect it is identical to Bernal- 
stacked graphite, but, as compared to Bernal-stacked 
graphite, there is a missing site at each end, and the 
missing sites are always free sites. We rename the sites 
Aj and Bj with aj and bj, so that aj and bj represents 
the chained site and the free site on the j-th layer, re- 
spectively. If ao = Bo, for example, the relation is 



Bj (j = even) 

Aj (j=odd) ' 

Aj (j — even) 

Bj (j = odd) ■ 



(G) 



The Hamiltonian is obtained by eliminating the miss- 
ing sites (bo and b^) in complete N + 1-layer Bernal 
graphite. As the system includes 2N atoms, there are 
2N eigenfunctions at a momentum p, each of which may 
be expressed as 



ipx 



(7) 



where x = (x, y) is the in-plane position and p = (p x ,p y ) 
is the in-plane momentum measured from K^. In this 
simple model with 70 and 71 , the energy bands are always 
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isotropic around p 
of p = 



0, and the dispersion is a function 



hi 



As we describe in the following, it is possible to classify 
the bases of the eigenfunctions into five categories: 

CI: chained, linear 
C2: chained, quadratic 
Fl: free, linear 
F2: free, quadratic 

F3: free, boundary-localized. (8) 

C and F represent chained and free sites, respectively. 1 
and 2 correspond to the linear and quadratic dispersion in 
the band structure. F3 gives a dispersion-less flat band in 
the isolated incomplete graphite, but it forms to a cubic 
(or higher order) band when the inter-section coupling is 
included. As an example, we present in Fig. 0] the band 
structure and the wave functions of incomplete graphite 
N = 8. 

The bases of CI and C2 are plane waves on the chained 
sites, which vanish at imaginary sites a_x and a^^i out- 
side the system. Explictly, this is defined as 



- sin qi (j 



N+2 



1) 



with quantized wave numbers 



N 



I = 1,2,- 



,N+l. 



(9) 



(10) 



The wavenumber qi = 7r/2, appearing only when N is 
even, is categorized in CI, and all the others in C2. 

The bases of Fl and F2 are plane waves on the free 
sites, but they have nodes at b\ and 6jv-i inside the sys- 
tem. They are written as 



o 



W=2 sin 9;'0' ~ x ) 



(11) 



with a different series of wave numbers 
, lir 



N -2' 



1=1,2,.-- ,JV-3, (12) 



where 0j is defined by e j — (p x ± ip y )/p when bj = Bj 
and Aj, respectively. Fl or F2 states do not exist when 
N < 3. The wavenumber q[ = ir/2 which appears when 
N is even (more than 4) is categorized in Fl, and all the 
others in F2. 

The group F3 is comprised of two states localized at 
&i or bjv-ii 



/l F3 (M 



hi) ' 
o 

5j.N-l 



(13) 



In the case N — 2, there is a single F3 state localized at 
the only free site. When N = 1, there are no F3 states. 

At p = 0, CI states and all F states (Fl, F2 and F3) 
are the exact eigenstates at zero-energy. Considering the 
perturbation in p for these degenerate states, the Hamil- 
tonian are approximately block diagonalized into blocks 
describing Cl+Fl, F2 and F3 states. CI and Fl are hy- 
bridized by a term linear in p, to form a monolayer-like 
Dirac cone, 



.Cl+Fl 



(P)»± 



N —2 
N + 2 



vp. 



(14) 



but with a reduced band velocity compared to mono- 
layer graphene (32j . Mixing with other states gives rise 
to a correction to Eq. (|14l) , but it is quadratic in p. The 
Cl+Fl band only appears when N is an even number, 
because otherwise CI or Fl state does not exist. The 
case of N — 2 is special, in that Fl does not exist and 
CI alone gives a flat band at zero energy. 
F2 states give low-energy quadratic bands, 



2 2 

_F2 / \ _ v P 



#(P)«- 



271 cos q[ 



(15) 



while F3 remains at zero energy independently of p, 



_F3 



(P) 



_F.3 



(p) = 0. 



(16) 



In Appendix [A] wc show that F2 and F3 states are actu- 
ally the approximate eigenstates with the eigen energies 
Eqs. (fT5|l and (fll)]) . respectively. F2 bands of q[ and 7r — q[ 
form a pair of electron and hole bands with the same band 
mass m* = 71 1 cosq[\/v 2 , analogous to the low-energy 
bands of bilayer graphene with the mass m* = 7 1 /(2i> 2 ) 
[3^ ]. Indeed, the 2x2 effective Hamiltonian spanned 
by these two F2 states is shown to be equivalent with 
the low-energy Hamiltonian for bilayer graphene by an 
appropriate unitary transformation. 

C2 states are the eigenstates on the chained sites at 
p = 0, with the non-zero eigen energy 271 cos qi. For 
p / 0, those states are hybridized with the free sites 
giving rise to quadratic dispersion for vp <?C 271 cos qi , 



-C2 



(p) ~ 27! cos qi + 



2 2 
v z p z 



271 cos qi 



1 



iV + 2 



sin q t 



(17) 



This band is electron-type and hole- type for < qi < it/2 
and 7r/2 < qi < n, respectively. 

In a section located at the end of the whole system as 
in Fig. OHc), one of missing sites, 60 or is restored, 
giving a different energy spectrum. When bo (bjsr) exists, 
the plane waves of free sites in Eq. (fTTj) have a node at 
6_i(6jv+i) (out of the system) instead of &i(6jv-i), so 
that the quantized wavenumbers on F states, Eq. 03 
are changed to 



In 



I = 1,2, 



,N- 



(18) 
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FIG. 4: (Left) Energy bands of incomplete Bernal graphite with N — 8. (Right) Some of wavefunctions from different 
categories at a small momentum along p x direction with Ip^l <C 71/11. White and black circles represent positive and negative 
wave amplitudes, respectively, and the size is the absolute magnitude. 



The dispersion of the F2 band in Eq. ([15]) changes ac- 
cording to the new set of q[. There are no changes for 
the C state wavenumber in Eq. (fTU)) . The linear band of 
Cl+Fl, Eq. (HI, is modified to 



£ 



Cl+Fl 



N 



N + 2 



Dp. 



(19) 



For F3 states, either (F3,L) or (F3,R), whichever is closer 
to the end of the whole system, does not exist any more. 

When both of 60 and bjy are restored, the partial sys- 
tem becomes a complete (N + l)-layer Bernal graphite, 
and qi and q[ become identical to Eq. ([TUl) . Then C2 
and F2 bands belonging to qi and ir — qi compose a sub- 
system equivalent to bilayer graphene, and the Cl+Fl 
band forms a linear band with band velocity equal to the 
monolayer's [H[27j]. 
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B. Coupling beyond the stacking faults 



FIG. 5: (a) Schematic of a stacking fault between the neigh- 
boring Bernal sections, (b) Structure of (Ni-i, 2, JVj_r). 



The wavefunctions of neighboring incomplete graphite 
sections are generally hybridized at p ^ via the bond 
between the ao site of one section and ajv of the other, as 
illustrated in Fig. [SJa). Particularly strong coupling oc- 
curs between the Cl+Fl band and the C2 band because 
they have significant wave amplitudes at the end sites oq 



and ajv, 



\fZ(*N) 



N + 2 



sin qi. 



(20) 
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The coupling is generally large at qi — ir/2 (CI state), 
and decreases as qi approaches or n. The F2 and F3 
bands, on the other hand, are mainly localized on free 
sites so that the mixing effect is generally small. In 
the following we consider details of the inter-section hy- 
bridization for each eigenstate class. 



1. Cl+Fl band 

If the sequence of integers {N\, iVjj, N3, ■ ■ ■ ,Nm) de- 
scribing the lengths JV» of each Bernal section con- 
tains any consecutive even numbers, then the Cl+Fl 
linear bands in those sections are mixed together 
through the coupling among CI states. Let us assume 
N s , Ng+i, • • • , N s+p are a series of consecutive even num- 
bers. We arrange the bases as \s, CI), |s,Fl), • • • , \s + 
p,Cl),\s + p, Fl), where \i, CI) represents the CI state 
of the section of iVj. The Hamiltonian matrix is written 
for this basis as 



H 



Cl+Fl _ 



fh s U s 

u\ h s +i u s+ i 

U l + 1 h s+2 U s+2 



,t 



(21) 



with 



f N l -2{l-5 i ) / vn] 

\ VTTi 



Ni + 2 



1 



U7r] 

(iV i + 2)(7V i+ i+2) V 



(22) 
(23) 



where 7r$ = 7r(7rT) when ao of the section s is A{B) site, 
and Si = 1 when the section Ni is located at the end of 
the entire system (i.e., i = 1 or M), and Si — otherwise. 

The Hamiltonian Eq. (|2"Tj) immediately leads to p + 1 
Dirac cones with generally different band velocities from 
the original. We found that the velocities are always 
equal to, or smaller than the monolayer's band veloc- 
ity v, and that velocity equal to v only appears when 
N s , ■ ■ ■ , N s+p covers the whole system, i.e., s = 1 and 
s + p = M. 

Fig. [5] shows an example of the band structure and 
the wavefunctions, calculated for the graphene stack 
(6,4,1,2). Dashed curves plotted for p < and E > 
represent the energy bands of the decomposed incom- 
plete graphites, N — 6,4,1, and 2. We see that the 
linear bands of the neighboring sections N = 6 and 4 are 
hybridized to give the 5th and 7th bands which have dif- 
ferent velocities from the original, while the linear band 
from N — 2 is kept almost intact (6th band). 



2. F2 band 

The coupling between F2 bands in the neighboring sec- 
tions is in proportional to p 3 , because the amplitude at 
ao and ajv is of the order of p, and the hopping between 
ao and ajv is vp. For small momenta vp <C 71, the cou- 
pling effect is negligible compared to the original F2 band 
energies oc p 2 , so that the wavefunction of F2 states are 
well localized in each single section. 

An exceptionally strong inter-section coupling occurs 
where a section of iVj = 2 exists in the middle of the 
system, as illustrated in Fig. [5](b). There CI states of 
Ni = 2 strongly hybridizes F2 states of the adjacent sec- 
tions i — 1 and i+X, and also two F3 states at 6jVi_i— 1 an d 
at b'l facing to the Ni = 2. The hybridized eigen energies 
are given by Eq. (|15[) with reconstructed wavenumber 
which are the solutions of 

sin [«rf(JV<_i - 1)] sin [q[{N l+1 - 1)] cos 2 q\ 

= cos fe'(AVi - 1)] cos [q[{N l+1 - 1)] sin 2 q\. (24) 

If the section i ± 1 is the end section, then N%±\ should 
be replaced with Ni±i + 2. The detailed derivation of 
Eqs. (f2"4")l is presented in the Appendix [Bl 



3. F3 band 

Similarly to F2, the coupling between F3 bands in 
neighboring sections is proportional to p 3 . Since the F3 
is originally a zero energy band in an isolated Bernal 
graphite, the coupling of 0(p 3 ) directly becomes the band 
dispersion when the inter-section coupling is introduced. 
Let us we consider F3 states \i, F3, R) and F3, L) in 
the neighboring sections i and which are facing each 
other across the stacking fault. We index the sites with 
a,j, bj and a'p b'j for the section i and i + 1, respectively, 
as illustrated in Fig.[5][a). The state \i, F3, R) is localized 
at fe/Vi— 1, and \i + 1,F3,L) is at b' x . When 6^-1 is an 
A site, the effective Hamiltonian in the basis of \i, F3, R) 
and \i + 1, F3, L) becomes, 



sn 



(eff) _ 











(25) 



leading to cubic dispersion of e = ±v 3 p 3 /j 2 . This ac- 
tually agrees with the cubic band of the stacking-fault 
bound states argued in Ref. [34[ . When &at ; _i is a B- 
site, 7r and 7T are interchanged in Eq. (|2"5")l . 

The Hamiltonian Eq. (I25|) is quite similar to that 
of low-energy states in rhombohedral trilayer graphene 
[24l l31 r 33] . If the stacking faults appear s times consec- 
utively, i.e., (• • • , Ni, 1, • • • ,1, N i+S , • • • ), then the states 
\i, F3, R) and \i + s, F3, L) forms a Hamiltonian equivalent 
to rhombohedral (s + 2)-layer, giving p s + 2 dispersion. In 
the multilayer (6,4,1,2), we have a pair of p 3 bands (2nd 
band) at the boundary between N = 6 and 4, and p A 
bands (1st band) between TV = 4 and 2, as shown in Fig. 

E 
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4. C2 band 

The intersection coupling modifies the band mass of 
C2 bands in Eq. ([T7| . The strength of coupling is given 
by the product of vp and the wave amplitudes at the end 
sites, Eq. (j2"0)) . The mass change is irrelevant when N is 
large (^> 1) or \qi\ is close to or it. When the energies 
of C2 bands in neighboring sections happen to coincide 
at p = 0, the two bands are strongly hybridized to give a 
linear p term in addition to the original quadratic term. 
This situation always occurs when the same qi is shared 
by the neighboring sections. In the example of Fig. HI we 
see that both of the sections N — 4 and 1 have a C2 state 
of qi = 7r/3, and those two bands are strongly hybridized 
to the 9th and 10th bands, which are linear at p = 0. 



IV. STATISTICAL PROPERTIES 
A. Linear band velocity and quadratic band mass 

The velocity of linear bands and the effective mass of 
quadratic bands are important parameters characteriz- 
ing the low-energy spectrum which can be probed exper- 
imentally, for example, by magneto-optical spectroscopy 
pol |4~0 | . In Bernal-rhombohedral mixed graphene multi- 
layer, there are several frequently appearing values of ve- 
locity and mass, each of which is associated with a single 
incomplete graphite section or particular combinations of 
them. 

For an isolated incomplete graphite, the velocity of lin- 
ear band (Cl+Fl) is given by y/(N- 2)/(N + 2)v for a 
middle section [Fig. El^b)] , and \fW/ (N + 2)v for an end 
section [Fig. 12c)]. The effective mass of the quadratic 
bands (F2) is \2('j 1 /v 2 ) cosg ; '|, where q\ is given by Eqs. 
(fl~2|) and (fl~8f for a middle section and an end section, re- 
spectively. Fig. Eta) lists the linear band velocity and the 
quadratic band mass in incomplete graphite sections with 
some different N's. Solid and dashed bars represent the 
values of middle sections and end sections, respectively. 

Fig. EJb) presents the same quantities of the mixed 
multilayer stacks of all possible configurations from sin- 
gle layer to 7 layers. There we derived the velocity and 
the mass from the numerically calculated eigenenergies 
of the original Hamiltonian Eq. (J3J) at small momenta. 
In multilayer (2,3), for example, we see that the linear 
band velocity is actually carried over from an isolated 
graphite of N = 2, and the quadratic band mass is from 
N = 3 (both of the end-section type) which are listed 
in Fig. [7ja). For the case where even numbers appear 
consecutively like (2,4), on the other hand, the linear- 
band velocities do not coincide with those of the isolated 
systems as argued in the previous section. The velocities 
of hybridized linear bands are then calculated by diag- 
onalizing Eq. (|21[) . Similarly, we see that the quadratic 
band mass does not match with isolated ones, when "2" 
appears somewhere in the middle, as (1,2,3). Then the 



(a) velocity 
(in units of v) 


configuration 


1/2 


(o,2,2,o)(o,4,6,o) 


i/Vs 


(o,4,o) (o,4,2,4,o) 


l/y/2 


(2,o)(o,6,o)(o,2,4,o)(o,2,2,2,o) 


a/375 


(0,8,0) 


v/273 


(4,o)(o,10,o)(o,4,6,o)(o,2,4,2,o)(o,4,2,4,o) 


V3/2 


(6,o)(2,2,o)(o,2,10,o)(o,6,8,o)(o,2,6,2,o) 




(b) mass (in 
units of 71 /2v 2 ) 


configuration 


(±l + V5)/2 


(5,.)(10,.)(.,7,.)(.,2,2,7,.)(.,5,2,5,.) 


1 


(3,.)(6,.)(.,5,.)(.,8,.)(.,3,2,3,.) 
(.,2,2,5,.)(.,3,2,6,.)(.2,2,8,.)(.,5,2,5,.) 




(4,.)(8,.)(1,2,1.)(.,6,.)(.,10,.)(.1,2,3,.) 
(.,1,2,7,.)(.,2,2,6,.)(.3,2,5,.)(.,4,2,4,.) 




(6,.)(.,8,.)(.1,2,4,.)(.,2,2,8,.)(.,4,2,7,.) 



TABLE I: (a) List of frequently-appearing velocity of the 
low-energy linear band. Symbol o in (o,a,b,o) represents an 
arbitrary sequence of numbers in which no even number comes 
next to a or b. (b) List of frequently-appearing effective mass 
of the low-energy quadratic band. Symbol • in (»,a,b,») rep- 
resents any sequences in which a number other than 2 comes 
next to a or b. 

hybridized mass is given by solving Eq. (1241) . 

The hybridized bands very often have the same band 
dispersions as a single incomplete graphite section. Ta- 
bles QJa) and b) list typical values of the linear band 
velocity and the quadratic band mass, respectively, with 
the corresponding configurations. In the velocity table, 
the sequence (o,a,6,o) represents the combination a, b ap- 
pearing in the middle of the multilayer, and (a,&,o) is one 
located at the end. The symbol o represents an arbitrary 
sequence of numbers in which no even number comes next 
to a or b (otherwise the linear bands are hybridized). In 
the mass table, similarly, • represents any sequences in 
which a number other than 2 comes next to a and b. 



B. Statistics of velocity and mass 

A realistic graphene multilayer system is not always 
a pure system with a single stacking arrangement, but 
is often composed of many small domains with different 
stacking structures [l(|. Nevertheless we expect to ob- 
serve strong signals from special band velocities or masses 
as argued above, because they are statistically likely to 
appear in a random configuration. Here we consider the 
distribution of band velocity and mass ensemble-averaged 
over possible Bernal-rhombohedral mixed configurations. 
Apart from tight-binding parameters 70 and 71 , the only 
parameter to characterize the system is p s , the probabil- 
ity to have an ABC stacking fault in every layer. Namely, 
the system is a pure Bernal multilayer when p s = 0, and 
pure rhombohedral multilayer when p s = 1. Experimen- 
tally, the value of p s is not well known, and it should 
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Graphene stack (6,4,1,2) 




vp/Yi 

FIG. 6: Band structure and some wavefunctions of the graphene stack (6, 4, 1, 2). Dashed curves plotted for p < and E > 
represent the energy bands of the decomposed incomplete graphites, N = 6, 4, 1, and 2 (iV = 6 and 2 are the end-section type). 

depend on the multilayer growing process. From the fact 
that rhombohedral structure is found by a few 10% out 
of natural graphite [7J, we expect p s is typically a few 
times 0.1. 

Figs. [5] and H] show the distribution of the linear band 
velocity and quadratic band mass in 100-layer systems 
averaged over 1000 different configurations, at p s = 0.1, 
0.2, and 0.4. The velocities and the masses are again 
obtained from the eigenenergies of the original Hamil- 
tonian Eq. §3§ at small momenta. The width of each 
histogram bin is 0.003w for velocity, and 0.006(7i/2i> 2 ) 
for mass. The height represents the average number of 
bands per a single 100-layer system. We actually observe 
the strong peaks at the frequently-appearing values listed 
in Table HI The relative heights between different peaks 
depends significantly on p s : In particular, the velocity of 
1/2 (in units of v) is almost absent in p s = 0.1, but at 
p s = 0.4 it is more than half of the tallest peak of 1/V2- 
It might be possible to estimate p s experimentally for a 
given graphite sample, by comparing the amplitudes of 
different peaks. 

We can approximately estimate the peak heights using 
simple probability calculation without taking the aver- 
age of the exact band structure. Let us consider a large 
multilayer stack composed of N tot layers, in which the 
stacking fault takes place with probability p s . The prob- 
ability for a given section (surrounded by stacking faults) 



to have length N is written as 

P(N)=p s (l-p s ) N - 1 . (26) 
The averaged length of a section is given by 

N = NP(N) = — , (27) 

N=l Ps 

and the average number of sections included in the whole 
system is 

M=-^=p s N tot . (28) 

Now we consider the expected number of a certain se- 
quence (o,a, o) existing in the whole system. When we 
look at a particular section, the probability that its length 
is a and the left- and right-neighboring sections are both 
odd numbers at the same time, is given by P ddP(a)Podd, 
where P odd = P(l) +P(3) +P(5) + • • • =l/(2-p s ). The 
number of the sequence (o, a, o) in the system is then ob- 
tained by multiplying this with the number of sections, 
M, namely, 

n(o, a, o) « Ps N tot P(a)P^ d . (29) 
For the end section (a, o), similarly, we have 

n(a,o)«2P(a)P odd , (30) 
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(a) Incompete graphite 

N 



velocity of linear band 
(in units of v) 

0.5 1 



mass of quadratic band 
(in units of Y./2V 2 ) 
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(b) General graphene stack 



f layers 



configuration 
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6 ABABABA 
1,5 ABCBCBC 
2,4 ABACACA 
3,3 ABABCBC 
1,1,4 ABCACAC 
1,2,3 ABCBABA 

1.3.2 ABCBCAC 

1.4.1 ABCBCBA 

2.1.3 ABACBCB 

2.2.2 ABACABA 
1,1,1,3 ABCABAB 
1,1,2,2 ABCACBC 

1.1.3.1 ABCACAB 

1.2.1.2 ABCBACA 

1.2.2.1 ABCBABC 

2.1.1.2 ABACBAB 
1,1,1,1,2 ABCABCB 
1,1,1,2,1 ABCABAC 
1,1,2,1,1 ABCACBA 

1,1,1,1,1,1 ABCABCA 
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FIG. 7: (a) Lists of the linear band velocity and the quadratic band mass in individual incomplete Bernal graphite sections 
from N = 1 to 14. Solid and dashed bars represent the values of middle sections and end sections, respectively, (b) List of the 
same quantities for multilayer stacks of all possible configurations from a single layer to 7 layers. 
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where the factor 2 comes from two ends of the whole 
system. For the combinate sequence (o,a, b, o), P(a) is 
just replaced with P(a)P(b) in above equations. For the 
case (»,a, •) appearing in the effective mass table, P oc id 
is replaced with 1 — P(2). 

The population of a given velocity/mass can be es- 
timated by summing up these numbers over all the 
corresponding sequences. For example, the number of 
bands with velocity l/\/3 is approximated by n(o, 4, o) + 
n(o, 4, 2,4, o), according to Table UJa). We calculate the 
populations of several representative velocities/masses in 
this manner, and show the estimated values as short hor- 
izontal bars in Figs. [8] and [9l While we took only a finite 
number of the configurations giving relatively large con- 
tributions, the approximation actually works quite well. 
We note that the peak diagrams like Figs. |S] and H] are 
not universal but depend on the total number of layers 
TVtot, because the contribution of the middle section is 
proportional to N tot as in Eq. (|2T))) , while that of the end 
section, Eq. (|50"|) . is not. In very large multilayer stacks 
with 7V tot >~ 1000, the end section component becomes 
negligible, and all the peaks just scale in proportional to 
AW 



C. Optical absorption spectrum 

Frequently appearing band structures give rise to strik- 
ing signals in the optical absorption spectrum. The op- 
tical absorption is related to the dynamical conductivity, 
which is given by within the linear response, 



/(£<*)- /M l<akl/3)| 2 



a,/3 



■ep 



, (31) 



where S is the area of the system, v x = &H/dp x is the 
velocity operator, S is the positive infinitesimal, f(e) is 
the Fermi distribution function, and \ot) and e a describe 
the eigenstate and the eigenenergy of the system. The 
transmission of light incident per pendicular to a two- 
dimensional system is given by [42J 



T = 



1 



2tt 



-a(u) 



(32) 



Here we calculate the dynamical conductivity of 100- 
layer graphites averaged over 1000 different configura- 
tions at fixed stacking- fault probability p s . We set the 
Fermi energy at the charge neutral point, ep = 0. Fig. 
ITOl shows the dynamical conductivity as a function of fre- 
quency uj calculated for p s = 0, 0.2 and 0.4. The peaks 
characterizing the spectrum mainly come from the transi- 
tion between e = and the band edge of the CI band, e = 
271 cos qi . The system at p s = is a pure Bernal-stacked 
100-layer graphite, and the absorption edge appears at 
hw = 271 cosg; with q t = Ztt/101 (I = 1,2, • • • , 100). In- 
deed, we see a regular series of small peaks for huj < 2ji . 
A major peak at hw — 71 comes from densely distributed 
absorption edges at qi ~ 0. 
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o 
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=> 
o 
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0.5 



Linear band velocity 
1 00 layers 

Ps = 0.1 1/V2 
1/V3 



/3/5 



—i \~ 

Ps = 0.2 



2/3 



I , J -lLi.I 



1/V3 



1/V2 



3/5 



1/2 



— 1 1— 

Ps = 0.4 



2/3 
s/l/2 



1/V3 1/V2 



1/2 



uL 



/3/5 



2/3 
v/3/2 



0.2 0.4 0.6 0.8 1 

velocity of linear band (in units of v) 

FIG. 8: Averaged distribution of the velocity of low- 
energy linear bands in 100-layer Bernal-rhombohedral mixed 
graphite with the stacking-fault probability p s = 0.1, 0.2, and 
0.4. Short horizontal bars attached to some major peaks are 
obtained by the approximation (see the text). 



For finite p s , on the other hand, the spectrum is 
rather regarded as a summation over isolated incomplete 
graphites of finite layers, as long as the total number 
of layers is much larger than the averaged length of the 
Bernal sections. The histogram in Fig. ITOl shows the dis- 
tribution of the energy of CI band edge e — 271 cos qi 
in a 100-layer system at p s = 0.2. We clearly see that 
the major peaks appearing in the dynamical conductivity 
correspond to the absorption edges of frequently appear- 
ing band structures. 



V. CONCLUSION 

We studied the electronic structures of general Bernal- 
rhombohedral mixed graphene multilayers with arbitrary 
configurations. We showed that every low-energy eigen- 
state is localized in a finite Bernal section bound by 
rhombohedral-type stacking faults, and the spectrum is 
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Quadratic band mass 
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Ps = 0.1 

(-1 + v/5)/2 
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Ps = 0.2 



(-1 + V5)/2 



1. 1 , . I j | I , r L. | l .i/.l^l.i. 



V'2 



(l + V§)/2 
1 



il t ,l | .JiJ.ilJ-L J]|l 



Ps = 0.4 



l 



(-l + V5)/2 



.1 1 i .T' 



(1 + VS)/2 
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0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 
mass of quadratic band (in units of y^/2\/ z ) 



1 00 layers 




p = 0.2 



15 



0.5 1 1.5 2 2.5 

Frequency fico & band edge energy 
(in units of y,) 

FIG. 10: Averaged dynamical conductivity of 100-layer 
Bernal-rhombohedral mixed graphites with stacking-fault 
probability p s = 0, 0.2 and 0.4. Histogram at the bottom 
is the averaged distribution of the band edges (band energies 
at zero momentum) in p a = 0.2. 
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FIG. 9: Averaged distribution of effective mass of low-energy 
quadratic bands similar to Fig. [8] 



well approximated by a collection of the spectra of in- 
dependent sections, categorized into linear, quadratic 
and cubic (or higher order) bands. We found that the 
ensemble-averaged spectrum is not smooth, but exhibits 
a number of discrete structures originating from finite 
Bernal sections or their combinations which are statis- 
tically likely to appear in a random configuration. In 
the low-energy region, in particular, there are frequently- 
appearing linear bands and quadratic bands with specific 
band velocities or curvatures. In the higher energy region 
away from the Dirac point, band edges are likely to ap- 
pear at some particular energies. Those discrete proper- 
ties may be detected in general graphite samples using ex- 
perimental techniques such as optical or magneto-optical 
absorption spectroscopy, and those observations would be 
useful to probe the stacking structure of graphene multi- 
layers. 



Appendix A: Eigen energies of F2 and F3 states 

Here we show that F2 and F3 states defined in the 
text are approximate eigen states with eigen energies Eqs. 
(|15l) and (|16p , respectively, in an independent incomplete 
Bernal graphite. Let us consider a TV-layer middle section 
shown in Fig. EJb), and assign the site indices a,j(J = 
0, 1, ■ • • , N) and bj(j = 1, • • ■ , N — 1) as in the figure. 
The Schrodinger equation on the layer j = 1, • • • , N is 
written as 

ef(b j )=vpe ie *f(a j ) (Al) 
sf( aj ) = vpe^fibj) + 7 i[/K--i) + /(ai+i)], (A2) 
and on j = and as 

e/(ao) = 7i/(ai), 

e/(ajv) = 7i/( a w-i), (A3) 

where e is the eigen energy. 

For the F2 state, the exact wave function is given by 
Eq. (fTTjl . plus a correction due to the first order pertur- 
bation in the momentum p. The wave amplitudes on the 
free sites are of 0th order and given as in Eq. (fTTj) , 



f(bj) = Ce"' sin ql(j - 1), 



(A4) 
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where C is a constant. The amplitudes on the chained 
sites are then derived from Eq. (|A1[) as 

f( aj ) = — Csmql(j-1). (A5) 
vp 

For j = 1, • • • , N — 1, we substitute above /(o 3 -) for Eq. 
to find 



e 2 = (271 cosg^e + v 2 p 2 



giving 



v 2 p 2 



2ji cos q\ 



(A6) 



(A7) 



as long as vp <C 271 cos q[ . For the end sites j = and 
j = N, on the other hand, Eq. (|A3|) yields to 



sin q\ = 0, 



(A8) 



which stands within 0(p 2 ) since e oc p 2 . 

For F3 states, it is straightforward to show that the 
exact eigenstates at finite p are explicitly written as 



/(&!) = !, /(« ) 



vpe 



7i 



-; for(F3,L) 



= 1, /(ajv) = - 



vpe 



7i 



-; for(F3,R) 
(A9) 



and the eigen energies are exactly zero. 



middle. Then the wave function is written as 



f(bj) =Ce ie *smq(j-l) 

f(b") = C"e* e " singCy-A^ + l) 



f(a' ) = = C, 



(Bl) 



where C, C',C" and g are quantities to be solved, and 
6j(6'j) is defined by e^'(e<) = (p x ±ip y )/p when 6^(6^) 
is i?-site and A-site, respectively. 

The amplitudes at chained sites a,j and a", which are 
linear in p, can be obtained in the exactly same way as 
in Appendix [X] as 

£ 

f{aj) = — Csin l(j ~ X ) 

f(a'!) = — C"sing(j - + 1), (B2) 
J vp 

where e is the eigen energy. The Schrodinger equations 
at aNi-x and Oq are given by 

vpe~* N *-if(a' ) +7i/( a Ni-i-i) = e/Ca^-i), 
«pe- l ^7(«' 2 ) + 7i/K) = e/(oo). (B3) 
By substituting the wave amplitudes, we have 



C = -a-jCfesing^i - 1) - 71 sing(^_ 1 - 2)] 
C" = — C" [esin (Z (^ i+1 - 1) - 7l sing(^ +1 - 2)] 



(B4) 



Appendix B: Inter-section coupling by Nt = 2 

We consider the hybridization of F2 and F3 states in 
a series of incomplete Bernal sections (JV»_i , 2, Afj+i), as 
illustrated in Fig.[5Jb), to derive the equation (|2"l)l for the 
reconstructed wavenumber q[. We assign the site indices 
(cii,bi), (a^b'j) and (a", 6") for the section iVj_i, 2, and 
N i+ i, respectively, as shown in Fig. [Sfb). We assume 
either JVf±a (out of the figure) is not 2, and then we can 
neglect the coupling with them. 

We consider a small momentum p <C Jx/v, and treat 
it as a perturbation. The 0th order wavefunction for a 
hybridized state is expressed as a linear combination of 

F2 and F3 states of the sections iV* 1 and iVj+i having a 

node at b\ and &jv +i -i' respectively, and CI state in the 



Also, the Schrodinger equations at a' and a' 2 , 

v V e i9 **-if{a Ni _,) +7i/(a'i) = e/K) 
«pe</(aj) + 71/K) = (B5) 



leads to 



C = —= [Csin g (7V j; _i - 1) + C"sxaq(N i+1 - 1)] .(B6) 
v2 

The wavenumber q and the relative wave amplitudes 
C'/C and C"/C are obtained by solving Eqs. (|B4l) and 
(|B6p . After some algebra, we obtain the equation for 
q = q'i as Eq. (|2"H) . While the argument above is only 
valid when iVj±i > 3, we can show that Eq. (|24p stand 
also when either or both Ni±\ is 1 or 2. 
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